"""

Created on Fri Mar 01 14:56:56 2013

Author: Josef Perktold
"""

import warnings

import numpy as np
from numpy.testing import (
    assert_allclose,
    assert_almost_equal,
    assert_array_less,
    assert_equal,
)
import pandas as pd
import pytest

import statsmodels.stats.proportion as smprop
from statsmodels.stats.proportion import (
    confint_proportions_2indep,
    multinomial_proportions_confint,
    power_proportions_2indep,
    proportion_confint,
    samplesize_proportions_2indep_onetail,
    score_test_proportions_2indep,
)
from statsmodels.stats.tests.results.results_proportion import (
    res_binom,
    res_binom_methods,
)
from statsmodels.tools.sm_exceptions import HypothesisTestWarning
from statsmodels.tools.testing import Holder

probci_methods = {
    "agresti_coull": "agresti-coull",
    "normal": "asymptotic",
    "beta": "exact",
    "wilson": "wilson",
    "jeffreys": "bayes",
}


@pytest.mark.parametrize("case", res_binom)
@pytest.mark.parametrize("method", probci_methods)
def test_confint_proportion(method, case):
    count, nobs = case
    idx = res_binom_methods.index(probci_methods[method])
    res_low = res_binom[case].ci_low[idx]
    res_upp = res_binom[case].ci_upp[idx]
    if np.isnan(res_low) or np.isnan(res_upp):
        pytest.skip("Skipping due to NaN value")
    if (count == 0 or count == nobs) and method == "jeffreys":
        # maybe a bug or different corner case definition
        pytest.skip("Skipping nobs 0 or count and jeffreys")
    if method == "jeffreys" and nobs == 30:
        # something is strange in extreme case e.g 0/30 or 1/30
        pytest.skip("Skipping nobs is 30 and jeffreys due to extreme case problem")
    ci = proportion_confint(count, nobs, alpha=0.05, method=method)
    # we impose that confint is in [0, 1]
    res_low = max(res_low, 0)
    res_upp = min(res_upp, 1)
    assert_almost_equal(ci, [res_low, res_upp], decimal=6, err_msg=repr(case) + method)


@pytest.mark.parametrize("method", probci_methods)
def test_confint_proportion_ndim(method):
    # check that it works with 1-D, 2-D and pandas

    count = np.arange(6).reshape(2, 3)
    nobs = 10 * np.ones((2, 3))

    count_pd = pd.DataFrame(count)
    nobs_pd = pd.DataFrame(nobs)

    ci_arr = proportion_confint(count, nobs, alpha=0.05, method=method)
    ci_pd = proportion_confint(count_pd, nobs_pd, alpha=0.05, method=method)
    assert_allclose(ci_arr, (ci_pd[0].values, ci_pd[1].values), rtol=1e-13)
    # spot checking one value
    ci12 = proportion_confint(count[1, 2], nobs[1, 2], alpha=0.05, method=method)
    assert_allclose((ci_pd[0].values[1, 2], ci_pd[1].values[1, 2]), ci12, rtol=1e-13)
    assert_allclose((ci_arr[0][1, 2], ci_arr[1][1, 2]), ci12, rtol=1e-13)

    # check that lists work as input
    ci_li = proportion_confint(count.tolist(), nobs.tolist(), alpha=0.05, method=method)
    assert_allclose(ci_arr, (ci_li[0], ci_li[1]), rtol=1e-13)

    # check pandas Series, 1-D
    ci_pds = proportion_confint(
        count_pd.iloc[0], nobs_pd.iloc[0], alpha=0.05, method=method
    )
    assert_allclose(
        (ci_pds[0].values, ci_pds[1].values),
        (ci_pd[0].values[0], ci_pd[1].values[0]),
        rtol=1e-13,
    )

    # check scalar nobs, verifying one value
    ci_arr2 = proportion_confint(count, nobs[1, 2], alpha=0.05, method=method)
    assert_allclose((ci_arr2[0][1, 2], ci_arr[1][1, 2]), ci12, rtol=1e-13)

    # check floating point values
    ci_arr2 = proportion_confint(count + 1e-4, nobs[1, 2], alpha=0.05, method=method)
    # should be close to values with integer values
    assert_allclose((ci_arr2[0][1, 2], ci_arr[1][1, 2]), ci12, rtol=1e-4)


def test_samplesize_confidenceinterval_prop():
    # consistency test for samplesize to achieve confidence_interval
    nobs = 20
    ci = smprop.proportion_confint(12, nobs, alpha=0.05, method="normal")
    res = smprop.samplesize_confint_proportion(12.0 / nobs, (ci[1] - ci[0]) / 2)
    assert_almost_equal(res, nobs, decimal=13)


def test_proportion_effect_size():
    # example from blog
    es = smprop.proportion_effectsize(0.5, 0.4)
    assert_almost_equal(es, 0.2013579207903309, decimal=13)


def test_confint_multinomial_proportions():
    from .results.results_multinomial_proportions import res_multinomial

    for (method, description), values in res_multinomial.items():
        cis = multinomial_proportions_confint(values.proportions, 0.05, method=method)
        assert_almost_equal(
            values.cis,
            cis,
            decimal=values.precision,
            err_msg=f'"{method}" method, {description}',
        )


def test_multinomial_proportions_errors():
    # Out-of-bounds values for alpha raise a ValueError
    for alpha in [-0.1, 0, 1, 1.1]:
        with pytest.raises(ValueError):
            multinomial_proportions_confint([5] * 50, alpha=alpha)

    with pytest.raises(ValueError):
        multinomial_proportions_confint(np.arange(50) - 1)
    # Any unknown method is reported.
    for method in ["unknown_method", "sisok_method", "unknown-glaz"]:
        with pytest.raises(NotImplementedError):
            (multinomial_proportions_confint([5] * 50, method=method))


def test_confint_multinomial_proportions_zeros():
    # test when a count is zero or close to zero
    # values from R MultinomialCI
    ci01 = np.array(
        [
            0.09364718,
            0.1898413,
            0.00000000,
            0.0483581,
            0.13667426,
            0.2328684,
            0.10124019,
            0.1974343,
            0.10883321,
            0.2050273,
            0.17210833,
            0.2683024,
            0.09870919,
            0.1949033,
        ]
    ).reshape(-1, 2)

    ci0 = np.array(
        [
            0.09620253,
            0.19238867,
            0.00000000,
            0.05061652,
            0.13924051,
            0.23542664,
            0.10379747,
            0.19998360,
            0.11139241,
            0.20757854,
            0.17468354,
            0.27086968,
            0.10126582,
            0.19745196,
        ]
    ).reshape(-1, 2)

    # the shifts are the differences between "LOWER(SG)"  "UPPER(SG)" and
    # "LOWER(C+1)" "UPPER(C+1)" in verbose printout
    # ci01_shift = np.array([0.002531008, -0.002515122])  # not needed
    ci0_shift = np.array([0.002531642, 0.002515247])

    p = [56, 0.1, 73, 59, 62, 87, 58]
    ci_01 = smprop.multinomial_proportions_confint(p, 0.05, method="sison_glaz")
    p = [56, 0, 73, 59, 62, 87, 58]
    ci_0 = smprop.multinomial_proportions_confint(p, 0.05, method="sison_glaz")

    assert_allclose(ci_01, ci01, atol=1e-5)
    assert_allclose(ci_0, np.maximum(ci0 - ci0_shift, 0), atol=1e-5)
    assert_allclose(ci_01, ci_0, atol=5e-4)


class CheckProportionMixin:
    def test_proptest(self):
        # equality of k-samples
        pt = smprop.proportions_chisquare(self.n_success, self.nobs, value=None)
        assert_almost_equal(pt[0], self.res_prop_test.statistic, decimal=13)
        assert_almost_equal(pt[1], self.res_prop_test.p_value, decimal=13)

        # several against value
        pt = smprop.proportions_chisquare(
            self.n_success, self.nobs, value=self.res_prop_test_val.null_value[0]
        )
        assert_almost_equal(pt[0], self.res_prop_test_val.statistic, decimal=13)
        assert_almost_equal(pt[1], self.res_prop_test_val.p_value, decimal=13)

        # one proportion against value
        pt = smprop.proportions_chisquare(
            self.n_success[0], self.nobs[0], value=self.res_prop_test_1.null_value
        )
        assert_almost_equal(pt[0], self.res_prop_test_1.statistic, decimal=13)
        assert_almost_equal(pt[1], self.res_prop_test_1.p_value, decimal=13)

    def test_pairwiseproptest(self):
        ppt = smprop.proportions_chisquare_allpairs(
            self.n_success, self.nobs, multitest_method=None
        )
        assert_almost_equal(ppt.pvals_raw, self.res_ppt_pvals_raw)
        ppt = smprop.proportions_chisquare_allpairs(
            self.n_success, self.nobs, multitest_method="h"
        )
        assert_almost_equal(ppt.pval_corrected(), self.res_ppt_pvals_holm)

        pptd = smprop.proportions_chisquare_pairscontrol(
            self.n_success, self.nobs, multitest_method="hommel"
        )
        assert_almost_equal(
            pptd.pvals_raw, ppt.pvals_raw[: len(self.nobs) - 1], decimal=13
        )

    def test_number_pairs_1493(self):
        ppt = smprop.proportions_chisquare_allpairs(
            self.n_success[:3], self.nobs[:3], multitest_method=None
        )

        assert_equal(len(ppt.pvals_raw), 3)
        idx = [0, 1, 3]
        assert_almost_equal(ppt.pvals_raw, self.res_ppt_pvals_raw[idx])


class TestProportion(CheckProportionMixin):
    def setup_method(self):
        self.n_success = np.array([73, 90, 114, 75])
        self.nobs = np.array([86, 93, 136, 82])

        self.res_ppt_pvals_raw = np.array(
            [
                0.00533824886503131,
                0.8327574849753566,
                0.1880573726722516,
                0.002026764254350234,
                0.1309487516334318,
                0.1076118730631731,
            ]
        )
        self.res_ppt_pvals_holm = np.array(
            [
                0.02669124432515654,
                0.8327574849753566,
                0.4304474922526926,
                0.0121605855261014,
                0.4304474922526926,
                0.4304474922526926,
            ]
        )

        res_prop_test = Holder()
        res_prop_test.statistic = 11.11938768628861
        res_prop_test.parameter = 3
        res_prop_test.p_value = 0.011097511366581344
        res_prop_test.estimate = np.array(
            [
                0.848837209302326,
                0.967741935483871,
                0.838235294117647,
                0.9146341463414634,
            ]
        ).reshape(4, 1, order="F")
        res_prop_test.null_value = """NULL"""
        res_prop_test.conf_int = """NULL"""
        res_prop_test.alternative = "two.sided"
        res_prop_test.method = (
            "4-sample test for equality of proportions without continuity correction"
        )
        res_prop_test.data_name = "smokers2 out of patients"
        self.res_prop_test = res_prop_test

        # > pt = prop.test(smokers2, patients, p=rep(c(0.9), 4), correct=FALSE)
        # > cat_items(pt, "res_prop_test_val.")
        res_prop_test_val = Holder()
        res_prop_test_val.statistic = np.array([13.20305530710751]).reshape(
            1, 1, order="F"
        )
        res_prop_test_val.parameter = np.array([4]).reshape(1, 1, order="F")
        res_prop_test_val.p_value = 0.010325090041836
        res_prop_test_val.estimate = np.array(
            [
                0.848837209302326,
                0.967741935483871,
                0.838235294117647,
                0.9146341463414634,
            ]
        ).reshape(4, 1, order="F")
        res_prop_test_val.null_value = np.array([0.9, 0.9, 0.9, 0.9]).reshape(
            4, 1, order="F"
        )
        res_prop_test_val.conf_int = """NULL"""
        res_prop_test_val.alternative = "two.sided"
        res_prop_test_val.method = (
            "4-sample test for given proportions without continuity correction"
        )
        res_prop_test_val.data_name = (
            "smokers2 out of patients, null probabilities rep(c(0.9), 4)"
        )
        self.res_prop_test_val = res_prop_test_val

        # > pt = prop.test(smokers2[1], patients[1], p=0.9, correct=FALSE)
        # > cat_items(pt, "res_prop_test_1.")
        res_prop_test_1 = Holder()
        res_prop_test_1.statistic = 2.501291989664086
        res_prop_test_1.parameter = 1
        res_prop_test_1.p_value = 0.113752943640092
        res_prop_test_1.estimate = 0.848837209302326
        res_prop_test_1.null_value = 0.9
        res_prop_test_1.conf_int = np.array([0.758364348004061, 0.9094787701686766])
        res_prop_test_1.alternative = "two.sided"
        res_prop_test_1.method = (
            "1-sample proportions test without continuity correction"
        )
        res_prop_test_1.data_name = (
            "smokers2[1] out of patients[1], null probability 0.9"
        )
        self.res_prop_test_1 = res_prop_test_1

    # GH 2969
    def test_default_values(self):
        count = np.array([5, 12])
        nobs = np.array([83, 99])
        stat, pval = smprop.proportions_ztest(count, nobs, value=None)
        assert_almost_equal(stat, -1.4078304151258787)
        assert_almost_equal(pval, 0.15918129181156992)

    # GH 2779
    def test_scalar(self):
        count = 5
        nobs = 83
        value = 0.05
        stat, pval = smprop.proportions_ztest(count, nobs, value=value)
        assert_almost_equal(stat, 0.392126026314)
        assert_almost_equal(pval, 0.694965098115)

        with pytest.raises(ValueError):
            smprop.proportions_ztest(count, nobs, value=None)


def test_binom_test():
    # > bt = binom.test(51,235,(1/6),alternative="less")
    # > cat_items(bt, "binom_test_less.")
    binom_test_less = Holder()
    binom_test_less.statistic = 51
    binom_test_less.parameter = 235
    binom_test_less.p_value = 0.982022657605858
    binom_test_less.conf_int = [0, 0.2659460862574313]
    binom_test_less.estimate = 0.2170212765957447
    binom_test_less.null_value = 1.0 / 6
    binom_test_less.alternative = "less"
    binom_test_less.method = "Exact binomial test"
    binom_test_less.data_name = "51 and 235"

    # > bt = binom.test(51,235,(1/6),alternative="greater")
    # > cat_items(bt, "binom_test_greater.")
    binom_test_greater = Holder()
    binom_test_greater.statistic = 51
    binom_test_greater.parameter = 235
    binom_test_greater.p_value = 0.02654424571169085
    binom_test_greater.conf_int = [0.1735252778065201, 1]
    binom_test_greater.estimate = 0.2170212765957447
    binom_test_greater.null_value = 1.0 / 6
    binom_test_greater.alternative = "greater"
    binom_test_greater.method = "Exact binomial test"
    binom_test_greater.data_name = "51 and 235"

    # > bt = binom.test(51,235,(1/6),alternative="t")
    # > cat_items(bt, "binom_test_2sided.")
    binom_test_2sided = Holder()
    binom_test_2sided.statistic = 51
    binom_test_2sided.parameter = 235
    binom_test_2sided.p_value = 0.0437479701823997
    binom_test_2sided.conf_int = [0.1660633298083073, 0.2752683640289254]
    binom_test_2sided.estimate = 0.2170212765957447
    binom_test_2sided.null_value = 1.0 / 6
    binom_test_2sided.alternative = "two.sided"
    binom_test_2sided.method = "Exact binomial test"
    binom_test_2sided.data_name = "51 and 235"

    alltests = [
        ("larger", binom_test_greater),
        ("smaller", binom_test_less),
        ("two-sided", binom_test_2sided),
    ]

    for alt, res0 in alltests:
        # only p-value is returned
        res = smprop.binom_test(51, 235, prop=1.0 / 6, alternative=alt)
        # assert_almost_equal(res[0], res0.statistic)
        assert_almost_equal(res, res0.p_value, decimal=13)

    # R binom_test returns Copper-Pearson confint
    ci_2s = smprop.proportion_confint(51, 235, alpha=0.05, method="beta")
    ci_low, ci_upp = smprop.proportion_confint(51, 235, alpha=0.1, method="beta")
    assert_almost_equal(ci_2s, binom_test_2sided.conf_int, decimal=13)
    assert_almost_equal(ci_upp, binom_test_less.conf_int[1], decimal=13)
    assert_almost_equal(ci_low, binom_test_greater.conf_int[0], decimal=13)


def test_binom_rejection_interval():
    # consistency check with binom_test
    # some code duplication but limit checks are different
    alpha = 0.05
    nobs = 200
    prop = 12.0 / 20
    alternative = "smaller"
    ci_low, ci_upp = smprop.binom_test_reject_interval(
        prop, nobs, alpha=alpha, alternative=alternative
    )
    assert_equal(ci_upp, nobs)
    pval = smprop.binom_test(ci_low, nobs, prop=prop, alternative=alternative)
    assert_array_less(pval, alpha)
    pval = smprop.binom_test(ci_low + 1, nobs, prop=prop, alternative=alternative)
    assert_array_less(alpha, pval)

    alternative = "larger"
    ci_low, ci_upp = smprop.binom_test_reject_interval(
        prop, nobs, alpha=alpha, alternative=alternative
    )
    assert_equal(ci_low, 0)
    pval = smprop.binom_test(ci_upp, nobs, prop=prop, alternative=alternative)
    assert_array_less(pval, alpha)
    pval = smprop.binom_test(ci_upp - 1, nobs, prop=prop, alternative=alternative)
    assert_array_less(alpha, pval)

    alternative = "two-sided"
    ci_low, ci_upp = smprop.binom_test_reject_interval(
        prop, nobs, alpha=alpha, alternative=alternative
    )
    pval = smprop.binom_test(ci_upp, nobs, prop=prop, alternative=alternative)
    assert_array_less(pval, alpha)
    pval = smprop.binom_test(ci_upp - 1, nobs, prop=prop, alternative=alternative)
    assert_array_less(alpha, pval)
    pval = smprop.binom_test(ci_upp, nobs, prop=prop, alternative=alternative)
    assert_array_less(pval, alpha)

    pval = smprop.binom_test(ci_upp - 1, nobs, prop=prop, alternative=alternative)
    assert_array_less(alpha, pval)


def test_binom_tost():
    # consistency check with two different implementation,
    # proportion_confint is tested against R
    # no reference case from other package available
    ci = smprop.proportion_confint(10, 20, method="beta", alpha=0.1)
    bt = smprop.binom_tost(10, 20, *ci)
    assert_almost_equal(bt, [0.05] * 3, decimal=12)

    ci = smprop.proportion_confint(5, 20, method="beta", alpha=0.1)
    bt = smprop.binom_tost(5, 20, *ci)
    assert_almost_equal(bt, [0.05] * 3, decimal=12)

    # vectorized, TODO: observed proportion = 0 returns nan
    ci = smprop.proportion_confint(np.arange(1, 20), 20, method="beta", alpha=0.05)
    bt = smprop.binom_tost(np.arange(1, 20), 20, ci[0], ci[1])
    bt = np.asarray(bt)
    assert_almost_equal(bt, 0.025 * np.ones(bt.shape), decimal=12)


def test_power_binom_tost():
    # comparison numbers from PASS manual
    p_alt = 0.6 + np.linspace(0, 0.09, 10)
    power = smprop.power_binom_tost(0.5, 0.7, 500, p_alt=p_alt, alpha=0.05)
    res_power = np.array(
        [0.9965, 0.9940, 0.9815, 0.9482, 0.8783, 0.7583, 0.5914, 0.4041, 0.2352, 0.1139]
    )
    assert_almost_equal(power, res_power, decimal=4)

    rej_int = smprop.binom_tost_reject_interval(0.5, 0.7, 500)
    res_rej_int = (269, 332)
    assert_equal(rej_int, res_rej_int)

    # TODO: actual alpha=0.0489  for all p_alt above

    # another case
    nobs = np.arange(20, 210, 20)
    power = smprop.power_binom_tost(0.4, 0.6, nobs, p_alt=0.5, alpha=0.05)
    res_power = np.array(
        [0.0, 0.0, 0.0, 0.0889, 0.2356, 0.3517, 0.4457, 0.6154, 0.6674, 0.7708]
    )
    # TODO: I currently do not impose power>=0, i.e np.maximum(power, 0)
    assert_almost_equal(np.maximum(power, 0), res_power, decimal=4)


def test_power_ztost_prop():
    power = smprop.power_ztost_prop(
        0.1, 0.9, 10, p_alt=0.6, alpha=0.05, discrete=True, dist="binom"
    )[0]
    assert_almost_equal(power, 0.8204, decimal=4)  # PASS example

    with warnings.catch_warnings():  # python >= 2.6
        warnings.simplefilter("ignore", HypothesisTestWarning)
        power = smprop.power_ztost_prop(
            0.4,
            0.6,
            np.arange(20, 210, 20),
            p_alt=0.5,
            alpha=0.05,
            discrete=False,
            dist="binom",
        )[0]

        res_power = np.array(
            [0.0, 0.0, 0.0, 0.0889, 0.2356, 0.4770, 0.5530, 0.6154, 0.7365, 0.7708]
        )
        # TODO: I currently do not impose power>=0, i.e np.maximum(power, 0)
        assert_almost_equal(np.maximum(power, 0), res_power, decimal=4)

        # with critval_continuity correction
        power = smprop.power_ztost_prop(
            0.4,
            0.6,
            np.arange(20, 210, 20),
            p_alt=0.5,
            alpha=0.05,
            discrete=False,
            dist="binom",
            variance_prop=None,
            continuity=2,
            critval_continuity=1,
        )[0]

        res_power = np.array(
            [0.0, 0.0, 0.0, 0.0889, 0.2356, 0.3517, 0.4457, 0.6154, 0.6674, 0.7708]
        )
        # TODO: I currently do not impose power>=0, i.e np.maximum(power, 0)
        assert_almost_equal(np.maximum(power, 0), res_power, decimal=4)

        power = smprop.power_ztost_prop(
            0.4,
            0.6,
            np.arange(20, 210, 20),
            p_alt=0.5,
            alpha=0.05,
            discrete=False,
            dist="binom",
            variance_prop=0.5,
            critval_continuity=1,
        )[0]

        res_power = np.array(
            [0.0, 0.0, 0.0, 0.0889, 0.2356, 0.3517, 0.4457, 0.6154, 0.6674, 0.7112]
        )
        # TODO: I currently do not impose power>=0, i.e np.maximum(power, 0)
        assert_almost_equal(np.maximum(power, 0), res_power, decimal=4)


def test_ztost():
    xfair = np.repeat([1, 0], [228, 762 - 228])

    # comparing to SAS last output at
    # http://support.sas.com/documentation/cdl/en/procstat/63104/HTML/default/viewer.htm#procstat_freq_sect028.htm
    # confidence interval for tost
    # generic ztost is moved to weightstats
    from statsmodels.stats.weightstats import zconfint, ztost

    ci01 = zconfint(xfair, alpha=0.1, ddof=0)
    assert_almost_equal(ci01, [0.2719, 0.3265], 4)
    res = ztost(xfair, 0.18, 0.38, ddof=0)

    assert_almost_equal(res[1][0], 7.1865, 4)
    assert_almost_equal(res[2][0], -4.8701, 4)
    assert_array_less(res[0], 0.0001)


def test_power_ztost_prop_norm():
    # regression test for normal distribution
    # from a rough comparison, the results and variations look reasonable
    with pytest.warns(HypothesisTestWarning):
        power = smprop.power_ztost_prop(
            0.4,
            0.6,
            np.arange(20, 210, 20),
            p_alt=0.5,
            alpha=0.05,
            discrete=False,
            dist="norm",
            variance_prop=0.5,
            continuity=0,
            critval_continuity=0,
        )[0]

    res_power = np.array(
        [
            0.0,
            0.0,
            0.0,
            0.11450013,
            0.27752006,
            0.41495922,
            0.52944621,
            0.62382638,
            0.70092914,
            0.76341806,
        ]
    )
    # TODO: I currently do not impose power>=0, i.e np.maximum(power, 0)
    assert_almost_equal(np.maximum(power, 0), res_power, decimal=4)

    # regression test for normal distribution
    with pytest.warns(HypothesisTestWarning):
        power = smprop.power_ztost_prop(
            0.4,
            0.6,
            np.arange(20, 210, 20),
            p_alt=0.5,
            alpha=0.05,
            discrete=False,
            dist="norm",
            variance_prop=0.5,
            continuity=1,
            critval_continuity=0,
        )[0]

    res_power = np.array(
        [
            0.0,
            0.0,
            0.02667562,
            0.20189793,
            0.35099606,
            0.47608598,
            0.57981118,
            0.66496683,
            0.73427591,
            0.79026127,
        ]
    )
    # TODO: I currently do not impose power>=0, i.e np.maximum(power, 0)
    assert_almost_equal(np.maximum(power, 0), res_power, decimal=4)

    # regression test for normal distribution
    with pytest.warns(HypothesisTestWarning):
        power = smprop.power_ztost_prop(
            0.4,
            0.6,
            np.arange(20, 210, 20),
            p_alt=0.5,
            alpha=0.05,
            discrete=True,
            dist="norm",
            variance_prop=0.5,
            continuity=1,
            critval_continuity=0,
        )[0]

    res_power = np.array(
        [
            0.0,
            0.0,
            0.0,
            0.08902071,
            0.23582284,
            0.35192313,
            0.55312718,
            0.61549537,
            0.66743625,
            0.77066806,
        ]
    )
    # TODO: I currently do not impose power>=0, i.e np.maximum(power, 0)
    assert_almost_equal(np.maximum(power, 0), res_power, decimal=4)

    # regression test for normal distribution
    with pytest.warns(HypothesisTestWarning):
        power = smprop.power_ztost_prop(
            0.4,
            0.6,
            np.arange(20, 210, 20),
            p_alt=0.5,
            alpha=0.05,
            discrete=True,
            dist="norm",
            variance_prop=0.5,
            continuity=1,
            critval_continuity=1,
        )[0]

    res_power = np.array(
        [
            0.0,
            0.0,
            0.0,
            0.08902071,
            0.23582284,
            0.35192313,
            0.44588687,
            0.61549537,
            0.66743625,
            0.71115563,
        ]
    )
    # TODO: I currently do not impose power>=0, i.e np.maximum(power, 0)
    assert_almost_equal(np.maximum(power, 0), res_power, decimal=4)

    # regression test for normal distribution
    with pytest.warns(HypothesisTestWarning):
        power = smprop.power_ztost_prop(
            0.4,
            0.6,
            np.arange(20, 210, 20),
            p_alt=0.5,
            alpha=0.05,
            discrete=True,
            dist="norm",
            variance_prop=None,
            continuity=0,
            critval_continuity=0,
        )[0]

    res_power = np.array(
        [
            0.0,
            0.0,
            0.0,
            0.0,
            0.15851942,
            0.41611758,
            0.5010377,
            0.5708047,
            0.70328247,
            0.74210096,
        ]
    )
    # TODO: I currently do not impose power>=0, i.e np.maximum(power, 0)
    assert_almost_equal(np.maximum(power, 0), res_power, decimal=4)


def test_proportion_ztests():
    # currently only consistency test with proportions chisquare
    # Note: alternative handling is generic

    res1 = smprop.proportions_ztest(15, 20.0, value=0.5, prop_var=0.5)
    res2 = smprop.proportions_chisquare(15, 20.0, value=0.5)
    assert_almost_equal(res1[1], res2[1], decimal=13)

    res1 = smprop.proportions_ztest(
        np.asarray([15, 10]), np.asarray([20.0, 20]), value=0, prop_var=None
    )
    res2 = smprop.proportions_chisquare(np.asarray([15, 10]), np.asarray([20.0, 20]))
    # test only p-value
    assert_almost_equal(res1[1], res2[1], decimal=13)

    # test with integers, issue #7603
    res1 = smprop.proportions_ztest(
        np.asarray([15, 10]), np.asarray([20, 50000]), value=0, prop_var=None
    )
    res2 = smprop.proportions_chisquare(np.asarray([15, 10]), np.asarray([20, 50000]))
    # test only p-value
    assert_almost_equal(res1[1], res2[1], decimal=13)
    assert_array_less(0, res2[-1][1])  # expected should be positive


def test_confint_2indep():
    # alpha = 0.05
    count1, nobs1 = 7, 34
    count2, nobs2 = 1, 34

    # result tables from Fagerland et al 2015
    """
    diff:
    Wald 0.029 0.32 0.29
    Agresti-Caffo 0.012 0.32 0.31
    Newcombe hybrid score 0.019 0.34 0.32
    Miettinen-Nurminen asymptotic score 0.028 0.34 0.31
    Santner-Snell exact unconditional -0.069 0.41 0.48
    Chan-Zhang exact unconditional 0.019 0.36 0.34
    Agresti-Min exact unconditional 0.024 0.35 0.33

    ratio:
    Katz log 0.91 54 4.08
    Adjusted log 0.92 27 3.38
    Inverse sinh 1.17 42 3.58
    Koopman asymptotic score 1.21 43 3.57
    Chan-Zhang 1.22 181 5.00
    Agresti-Min 1.15 89 4.35

    odds-ratio
    Woolf logit 0.99 74 4.31
    Gart adjusted logit 0.98 38 3.65
    Independence-smoothed logit 0.99 60 4.11
    Cornfield exact conditional 0.97 397 6.01
    Cornfield mid-p 1.19 200 5.12
    Baptista-Pike exact conditional 1.00 195 5.28
    Baptista-Pike mid-p 1.33 99 4.31
    Agresti-Min exact unconditional 1.19 72 4.10
    """  # pylint: disable=W0105
    ci = confint_proportions_2indep(
        count1, nobs1, count2, nobs2, method="newcomb", compare="diff", alpha=0.05
    )
    # one decimal to upp added from regression result
    assert_allclose(ci, [0.019, 0.340], atol=0.005)
    ci = confint_proportions_2indep(
        count1, nobs1, count2, nobs2, method="wald", compare="diff", alpha=0.05
    )
    assert_allclose(ci, [0.029, 0.324], atol=0.005)
    ci = confint_proportions_2indep(
        count1, nobs1, count2, nobs2, method="agresti-caffo", compare="diff", alpha=0.05
    )
    assert_allclose(ci, [0.012, 0.322], atol=0.005)
    ci = confint_proportions_2indep(
        count1, nobs1, count2, nobs2, compare="diff", method="score", correction=True
    )
    assert_allclose(ci, [0.028, 0.343], rtol=0.03)

    # ratio
    ci = confint_proportions_2indep(
        count1, nobs1, count2, nobs2, compare="ratio", method="log"
    )
    assert_allclose(ci, [0.91, 54], rtol=0.01)
    ci = confint_proportions_2indep(
        count1, nobs1, count2, nobs2, compare="ratio", method="log-adjusted"
    )
    assert_allclose(ci, [0.92, 27], rtol=0.01)
    ci = confint_proportions_2indep(
        count1, nobs1, count2, nobs2, compare="ratio", method="score", correction=False
    )
    assert_allclose(ci, [1.21, 43], rtol=0.01)

    # odds-ratio
    ci = confint_proportions_2indep(
        count1, nobs1, count2, nobs2, compare="or", method="logit"
    )
    assert_allclose(ci, [0.99, 74], rtol=0.01)
    ci = confint_proportions_2indep(
        count1, nobs1, count2, nobs2, compare="or", method="logit-adjusted"
    )
    assert_allclose(ci, [0.98, 38], rtol=0.01)
    ci = confint_proportions_2indep(
        count1, nobs1, count2, nobs2, compare="or", method="logit-smoothed"
    )
    assert_allclose(ci, [0.99, 60], rtol=0.01)
    ci = confint_proportions_2indep(
        count1,
        nobs1,
        count2,
        nobs2,
        compare="odds-ratio",
        method="score",
        correction=True,
    )
    # regression test
    assert_allclose(ci, [1.246622, 56.461576], rtol=0.01)


def test_confint_2indep_propcis():
    # unit tests compared to R package PropCis
    # alpha = 0.05
    count1, nobs1 = 7, 34
    count2, nobs2 = 1, 34

    # > library(PropCIs)
    # > diffscoreci(7, 34, 1, 34, 0.95)
    ci = 0.0270416, 0.3452912
    ci1 = confint_proportions_2indep(
        count1, nobs1, count2, nobs2, compare="diff", method="score", correction=True
    )
    assert_allclose(ci1, ci, atol=0.002)  # lower agreement (iterative)
    # > wald2ci(7, 34, 1, 34, 0.95, adjust="AC")
    ci = 0.01161167, 0.32172166
    ci1 = confint_proportions_2indep(
        count1, nobs1, count2, nobs2, compare="diff", method="agresti-caffo"
    )
    assert_allclose(ci1, ci, atol=6e-7)
    # > wald2ci(7, 34, 1, 34, 0.95, adjust="Wald")
    ci = 0.02916942, 0.32377176
    ci1 = confint_proportions_2indep(
        count1, nobs1, count2, nobs2, compare="diff", method="wald", correction=False
    )
    assert_allclose(ci1, ci, atol=6e-7)

    # > orscoreci(7, 34, 1, 34, 0.95)
    ci = 1.246309, 56.486130
    ci1 = confint_proportions_2indep(
        count1,
        nobs1,
        count2,
        nobs2,
        compare="odds-ratio",
        method="score",
        correction=True,
    )
    assert_allclose(ci1, ci, rtol=5e-4)  # lower agreement (iterative)

    # > riskscoreci(7, 34, 1, 34, 0.95)
    ci = 1.220853, 42.575718
    ci1 = confint_proportions_2indep(
        count1, nobs1, count2, nobs2, compare="ratio", method="score", correction=False
    )
    assert_allclose(ci1, ci, atol=6e-7)


def test_score_test_2indep():
    # this does not verify the statistic and pvalue yet
    count1, nobs1 = 7, 34
    count2, nobs2 = 1, 34

    for co in ["diff", "ratio", "or"]:
        res = score_test_proportions_2indep(count1, nobs1, count2, nobs2, compare=co)
        assert_allclose(res.prop1_null, res.prop2_null, rtol=1e-10)

        # check that equality case is handled
        val = 0 if co == "diff" else 1.0
        s0, pv0 = score_test_proportions_2indep(
            count1, nobs1, count2, nobs2, compare=co, value=val, return_results=False
        )[:2]
        s1, pv1 = score_test_proportions_2indep(
            count1,
            nobs1,
            count2,
            nobs2,
            compare=co,
            value=val + 1e-10,
            return_results=False,
        )[:2]
        assert_allclose(s0, s1, rtol=1e-8)
        assert_allclose(pv0, pv1, rtol=1e-8)
        s1, pv1 = score_test_proportions_2indep(
            count1,
            nobs1,
            count2,
            nobs2,
            compare=co,
            value=val - 1e-10,
            return_results=False,
        )[:2]
        assert_allclose(s0, s1, rtol=1e-8)
        assert_allclose(pv0, pv1, rtol=1e-8)


def test_test_2indep():
    # this checks the pvalue of the hypothesis test at value equal to the
    # confidence limit
    alpha = 0.05
    count1, nobs1 = 7, 34
    count2, nobs2 = 1, 34

    methods_both = [
        ("diff", "agresti-caffo"),
        # ('diff', 'newcomb'),  # only confint
        ("diff", "score"),
        ("diff", "wald"),
        ("ratio", "log"),
        ("ratio", "log-adjusted"),
        ("ratio", "score"),
        ("odds-ratio", "logit"),
        ("odds-ratio", "logit-adjusted"),
        ("odds-ratio", "logit-smoothed"),
        ("odds-ratio", "score"),
    ]

    for co, method in methods_both:
        low, upp = confint_proportions_2indep(
            count1,
            nobs1,
            count2,
            nobs2,
            compare=co,
            method=method,
            alpha=alpha,
            correction=False,
        )

        res = smprop.test_proportions_2indep(
            count1,
            nobs1,
            count2,
            nobs2,
            value=low,
            compare=co,
            method=method,
            correction=False,
        )
        assert_allclose(res.pvalue, alpha, atol=1e-10)

        res = smprop.test_proportions_2indep(
            count1,
            nobs1,
            count2,
            nobs2,
            value=upp,
            compare=co,
            method=method,
            correction=False,
        )
        assert_allclose(res.pvalue, alpha, atol=1e-10)

        _, pv = smprop.test_proportions_2indep(
            count1,
            nobs1,
            count2,
            nobs2,
            value=upp,
            compare=co,
            method=method,
            alternative="smaller",
            correction=False,
            return_results=False,
        )
        assert_allclose(pv, alpha / 2, atol=1e-10)

        _, pv = smprop.test_proportions_2indep(
            count1,
            nobs1,
            count2,
            nobs2,
            value=low,
            compare=co,
            method=method,
            alternative="larger",
            correction=False,
            return_results=False,
        )
        assert_allclose(pv, alpha / 2, atol=1e-10)

    # test Miettinen/Nurminen small sample correction
    co, method = "ratio", "score"
    low, upp = confint_proportions_2indep(
        count1,
        nobs1,
        count2,
        nobs2,
        compare=co,
        method=method,
        alpha=alpha,
        correction=True,
    )

    res = smprop.test_proportions_2indep(
        count1,
        nobs1,
        count2,
        nobs2,
        value=low,
        compare=co,
        method=method,
        correction=True,
    )
    assert_allclose(res.pvalue, alpha, atol=1e-10)


def test_equivalence_2indep():
    # this checks the pvalue of the equivalence test at value equal to the
    # confidence limit
    alpha = 0.05
    count1, nobs1 = 7, 34
    count2, nobs2 = 1, 34
    count1v, nobs1v = [7, 1], 34
    count2v, nobs2v = [1, 7], 34

    methods_both = [
        ("diff", "agresti-caffo"),
        # ('diff', 'newcomb'),  # only confint
        ("diff", "score"),
        ("diff", "wald"),
        ("ratio", "log"),
        ("ratio", "log-adjusted"),
        ("ratio", "score"),
        ("odds-ratio", "logit"),
        ("odds-ratio", "logit-adjusted"),
        ("odds-ratio", "logit-smoothed"),
        ("odds-ratio", "score"),
    ]

    for co, method in methods_both:
        low, upp = confint_proportions_2indep(
            count1,
            nobs1,
            count2,
            nobs2,
            compare=co,
            method=method,
            alpha=2 * alpha,
            correction=False,
        )

        # Note: test should have only one margin at confint
        res = smprop.tost_proportions_2indep(
            count1,
            nobs1,
            count2,
            nobs2,
            low,
            upp * 1.05,
            compare=co,
            method=method,
            correction=False,
        )
        assert_allclose(res.pvalue, alpha, atol=1e-10)

        res = smprop.tost_proportions_2indep(
            count1,
            nobs1,
            count2,
            nobs2,
            low * 0.95,
            upp,
            compare=co,
            method=method,
            correction=False,
        )
        assert_allclose(res.pvalue, alpha, atol=1e-10)

        # vectorized
        if method == "logit-smoothed":
            # not correctly vectorized
            return
        res = smprop.tost_proportions_2indep(
            count1v,
            nobs1v,
            count2v,
            nobs2v,
            low * 0.95,
            upp,
            compare=co,
            method=method,
            correction=False,
        )
        assert_allclose(res.pvalue[0], alpha, atol=1e-10)


def test_score_confint_koopman_nam():

    # example Koopman, based on Nam 1995

    x0, n0 = 16, 80
    x1, n1 = 36, 40
    # x = x0 + x1
    # n = n0 + n1
    # p0 = x0 / n0
    # p1 = x1 / n1

    results_nam = Holder()
    results_nam.p0_roots = [0.1278, 0.2939, 0.4876]
    results_nam.conf_int = [2.940, 7.152]

    res = smprop._confint_riskratio_koopman(x1, n1, x0, n0, alpha=0.05)

    assert_allclose(res._p_roots, results_nam.p0_roots, atol=4)
    assert_allclose(res.confint, results_nam.conf_int, atol=3)

    table = [67, 9, 7, 16]  # [67, 7, 9, 16]
    resp = smprop._confint_riskratio_paired_nam(table, alpha=0.05)
    # TODO: currently regression test, need verified results
    ci_old = [0.917832, 1.154177]
    assert_allclose(resp.confint, ci_old, atol=3)


def test_power_2indep():
    # test against R
    pow_ = power_proportions_2indep(-0.25, 0.75, 76.70692)
    assert_allclose(pow_.power, 0.9, atol=1e-8)

    n = samplesize_proportions_2indep_onetail(
        -0.25, 0.75, 0.9, ratio=1, alpha=0.05, value=0, alternative="two-sided"
    )
    assert_allclose(n, 76.70692, atol=1e-5)

    power_proportions_2indep(-0.25, 0.75, 62.33551, alternative="smaller")
    assert_allclose(pow_.power, 0.9, atol=1e-8)

    pow_ = power_proportions_2indep(0.25, 0.5, 62.33551, alternative="smaller")
    assert_array_less(pow_.power, 0.05)

    pow_ = power_proportions_2indep(
        0.25, 0.5, 62.33551, alternative="larger", return_results=False
    )
    assert_allclose(pow_, 0.9, atol=1e-8)

    pow_ = power_proportions_2indep(-0.15, 0.65, 83.4373, return_results=False)
    assert_allclose(pow_, 0.5, atol=1e-8)

    n = samplesize_proportions_2indep_onetail(
        -0.15, 0.65, 0.5, ratio=1, alpha=0.05, value=0, alternative="two-sided"
    )

    assert_allclose(n, 83.4373, atol=0.05)

    # Stata example
    from statsmodels.stats.power import normal_sample_size_one_tail

    res = power_proportions_2indep(-0.014, 0.015, 550, ratio=1.0)
    assert_allclose(res.power, 0.7415600, atol=1e-7)
    n = normal_sample_size_one_tail(
        -0.014, 0.7415600, 0.05 / 2, std_null=res.std_null, std_alternative=res.std_alt
    )
    assert_allclose(n, 550, atol=0.05)
    n2 = samplesize_proportions_2indep_onetail(
        -0.014, 0.015, 0.7415600, ratio=1, alpha=0.05, value=0, alternative="two-sided"
    )
    assert_allclose(n2, n, rtol=1e-13)

    # with nobs ratio != 1
    # note Stata has reversed ratio compared to ours, see #8049
    pwr_st = 0.7995659211532175
    n = 154
    res = power_proportions_2indep(-0.1, 0.2, n, ratio=2.0)
    assert_allclose(res.power, pwr_st, atol=1e-7)

    n2 = samplesize_proportions_2indep_onetail(-0.1, 0.2, pwr_st, ratio=2)
    assert_allclose(n2, n, rtol=1e-4)


@pytest.mark.parametrize("count", np.arange(10, 90, 5))
@pytest.mark.parametrize("method", list(probci_methods.keys()) + ["binom_test"])
def test_ci_symmetry(count, method):
    n = 100
    a = proportion_confint(count, n, method=method)
    b = proportion_confint(n - count, n, method=method)
    assert_allclose(np.array(a), 1.0 - np.array(b[::-1]))


@pytest.mark.parametrize("nobs", [47, 50])
@pytest.mark.parametrize("count", np.arange(48))
@pytest.mark.parametrize("array_like", [False, True])
def test_ci_symmetry_binom_test(nobs, count, array_like):
    _count = [count] * 3 if array_like else count
    nobs_m_count = [nobs - count] * 3 if array_like else nobs - count
    a = proportion_confint(_count, nobs, method="binom_test")
    b = proportion_confint(nobs_m_count, nobs, method="binom_test")
    assert_allclose(np.array(a), 1.0 - np.array(b[::-1]))


def test_int_check():
    # integer values are required only if method="binom_test"
    with pytest.raises(ValueError):
        proportion_confint(10.5, 20, method="binom_test")
    with pytest.raises(ValueError):
        proportion_confint(10, 20.5, method="binom_test")
    with pytest.raises(ValueError):
        proportion_confint(np.array([10.3]), 20, method="binom_test")

    a = proportion_confint(21.0, 47, method="binom_test")
    b = proportion_confint(21, 47, method="binom_test")
    c = proportion_confint(21, 47.0, method="binom_test")
    assert_allclose(a, b)
    assert_allclose(a, c)


@pytest.mark.parametrize("count", np.arange(10, 90, 5))
@pytest.mark.parametrize("method", list(probci_methods.keys()) + ["binom_test"])
def test_ci_symmetry_array(count, method):
    n = 100
    a = proportion_confint([count, count], n, method=method)
    b = proportion_confint([n - count, n - count], n, method=method)
    assert_allclose(np.array(a), 1.0 - np.array(b[::-1]))


@pytest.mark.parametrize(
    "count, nobs, alternative, expected",
    [
        (1, 1, "two-sided", [1.0, 1.0]),
        (0, 1, "two-sided", [0.0, 0.0]),
        (1, 1, "larger", [0.0, 1.0]),
        (0, 1, "larger", [0.0, 0.0]),
        (1, 1, "smaller", [1.0, 1.0]),
        (0, 1, "smaller", [0.0, 1.0]),
    ],
)
def test_proportion_confint_edge(count, nobs, alternative, expected):
    # test against R 4.3.3 / DescTools / BinomCI and binom.test
    ci_low, ci_upp = proportion_confint(count, nobs, alternative=alternative)
    assert_allclose(np.array([ci_low, ci_upp]), np.array(expected))


@pytest.mark.parametrize("count", [100])
@pytest.mark.parametrize("nobs", [200])
@pytest.mark.parametrize(
    "method, alpha, alternative, expected",
    [
        # two-sided
        ("normal", 0.05, "two-sided", [0.4307048, 0.5692952]),
        ("normal", 0.01, "two-sided", [0.4089307, 0.5910693]),
        ("normal", 0.10, "two-sided", [0.4418456, 0.5581544]),
        ("agresti_coull", 0.05, "two-sided", [0.4313609, 0.5686391]),
        ("beta", 0.05, "two-sided", [0.4286584, 0.5713416]),
        ("wilson", 0.05, "two-sided", [0.4313609, 0.5686391]),
        ("jeffreys", 0.05, "two-sided", [0.4311217, 0.5688783]),
        # larger
        ("normal", 0.05, "larger", [0.0, 0.5581544]),
        ("agresti_coull", 0.05, "larger", [0.0, 0.557765]),
        ("beta", 0.05, "larger", [0.0, 0.560359]),
        ("wilson", 0.05, "larger", [0.0, 0.557765]),
        ("jeffreys", 0.05, "larger", [0.0, 0.5578862]),
        # smaller
        ("normal", 0.05, "smaller", [0.4418456, 1.0]),
        ("agresti_coull", 0.05, "smaller", [0.442235, 1.0]),
        ("beta", 0.05, "smaller", [0.439641, 1.0]),
        ("wilson", 0.05, "smaller", [0.442235, 1.0]),
        ("jeffreys", 0.05, "smaller", [0.4421138, 1.0]),
    ],
)
def test_proportion_confint(count, nobs, method, alpha, alternative, expected):
    # test against R 4.3.3 / DescTools / BinomCI and binom.test
    ci_low, ci_upp = proportion_confint(
        count, nobs, alpha=alpha, method=method, alternative=alternative
    )
    assert_allclose(np.array([ci_low, ci_upp]), np.array(expected))


@pytest.mark.parametrize("count", [100])
@pytest.mark.parametrize("nobs", [200])
@pytest.mark.parametrize(
    "alpha, alternative, expected, rtol",
    [
        # binom.exact with method="two.sided" returns values with a precision of up to 4 decimal places.
        (0.01, "two-sided", [0.4094, 0.5906], 1e-4),
        (0.05, "two-sided", [0.4299, 0.5701], 1e-4),
        (0.10, "two-sided", [0.44, 0.56], 1e-4),
        (0.01, "larger", [0.0, 0.5840449], 1e-7),
        (0.05, "larger", [0.0, 0.560359], 1e-7),
        (0.10, "larger", [0.0, 0.5476422], 1e-7),
        (0.01, "smaller", [0.4159551, 1.0], 1e-7),
        (0.05, "smaller", [0.439641, 1.0], 1e-7),
        (0.10, "smaller", [0.4523578, 1.0], 1e-7),
    ],
)
def test_proportion_confint_binom_test(count, nobs, alpha, alternative, expected, rtol):
    # test against R 4.3.3 / DescTools / exactci
    ci_low, ci_upp = proportion_confint(
        count, nobs, alpha=alpha, method="binom_test", alternative=alternative
    )
    assert_allclose(np.array([ci_low, ci_upp]), np.array(expected), rtol=rtol)
